Exemple modèle d’occupancy dynamique
Olivier Gimenez
19/10/2020
1. Lecture des données
On charge le package unmarked pour ajuster les modèles dynamique d’occupancy.
library(unmarked)Les données sont celles du stage de Valentin, jusqu’à 2017 si je ne m’abuse, mises à jour entre Valentin, Julie et Christophe, en particulier sur l’effort. On a les données de détections/non-détections dans l’objet y, de dimensions le nombre de cellules de la grille utilisée, le nombre d’occasions secondaires (les mois d’hiver) et le nombre d’années. A noter que pur prendre en compte l’effort de prospection nul pour certaines années sur certaines cellules, on a dans les obs des NA.
load("dat/WS.RData")
dim(y)## [1] 3547 4 232. Formatage des données
On met les données de détection/non-détection au format unmarked avec occasions primaire et secondaires.
Yter <- y[,,1]
for (i in 2:dim(y)[3]){
Yter <- cbind(Yter,y[,,i])
}
dim(Yter) # 3547 x 92 ou 92 est 4 occ x 23 ans## [1] 3547 92On vérifie que tous les sites ont au moins 1 obs, et on supprime les cellules sur lesquelles aucun échantillonnage n’a jamais été fait dans les données. On fait pareil pour les covariables.
ind <- apply(Yter, 1, function(x) all(is.na(x)))
sum(ind)## [1] 97Yter <- Yter[!ind, ]
y <- y[!ind,, ]
mask <- mask[!ind,]
dim(CovEnvbis) # covariables environnementales## [1] 3547 7CovEnvbis <- CovEnvbis[!ind, ]
dim(Sbis) # autocorrélation spatiale à courte distance = 10km = 8 voisins proches ## [1] 3547 23Sbis <- Sbis[!ind,]
dim(Lbis) # autocorrélation spatiale à longue distance = 150 km## [1] 3547 23Lbis <- Lbis[!ind,]
dim(PO) # la fameuse pression d'observation ## [1] 3547 23PO <- PO[!ind,]On rassemble les covariables dont les valeurs diffèrent selon la cellule uniquement.
sites.covs <- data.frame(
forest = CovEnvbis[,"p_forest"],
agr = CovEnvbis[,"p_agri"],
rock = CovEnvbis[,"p_rock"],
halt = CovEnvbis[,"p_halt"],
alt = CovEnvbis[,"p_alti"],
dbarr = CovEnvbis[,"p_dbarr"],
acc = CovEnvbis[,"p_road"])On rassemble les covariables dont les valeurs diffèrent selon la cellule et l’année.
year <- matrix(rep(c('01',
'02',
'03',
'04',
'05',
'06',
'07',
'08',
'09',
'10',
'11',
'12',
'13',
'14',
'15',
'16',
'17',
'18',
'19',
'20',
'21',
'22',
'23'), nrow(Yter)),
nrow = nrow(Yter),
ncol = 23,
byrow = TRUE)
trendyear <- matrix(rep(1:23,nrow(Yter)),
nrow=nrow(Yter),
ncol=23,
byrow=TRUE)
yearly.site.covs <- list(
PO = PO[,1:23], # effort echantillonnage (prospection theorique)
year = year, # effet annee
trendyear = trendyear, # tendance lineaire annee
SDAC = Sbis[,1:23], # autocorrelation courte distance
LDAC = Lbis[,1:23]) # autocorrelation longue distanceOn rassemble les covariables dont les valeurs diffèrent selon les occasions secondaires.
occ <- array(0,dim(y))
for (i in 1:dim(y)[1]){
for (j in 1:dim(y)[3]){
occ[i,1:dim(y)[2],j] <- c('1','2','3','4') # dec, jan, fev, mar
if (mask[i,j]) occ[i,1:dim(y)[2],j] <- NA
}
}
occbis <- occ[,,1]
for (i in 2:dim(y)[3]){
occbis <- cbind(occbis,occ[,,i])
}
dim(occbis) # 3450 x 92 ou 92 est 4 occ x 23 ans## [1] 3450 92obs.covs <- list(OCC = occbis)On crée un data.frame (un objet fourre-tout dans R) pour passer dans le package unmarked.
umf <- unmarkedMultFrame(y = Yter,
siteCovs = sites.covs,
yearlySiteCovs = yearly.site.covs,
obsCovs = obs.covs,
numPrimary = 23)En résumé, on a quoi dans ces données?
summary(umf)## unmarkedFrame Object
##
## 3450 sites
## Maximum number of observations per site: 92
## Mean number of observations per site: 57.91
## Number of primary survey periods: 23
## Number of secondary survey periods: 4
## Sites with at least one detection: 547
##
## Tabulation of y observations:
## 0 1 2 <NA>
## 196086 3215 503 117596
##
## Site-level covariates:
## forest agr rock halt
## Min. :-1.218838 Min. :-1.50457 Min. :-0.253118 Min. :-0.172297
## 1st Qu.:-0.939161 1st Qu.:-0.83057 1st Qu.:-0.253118 1st Qu.:-0.172297
## Median :-0.102772 Median : 0.04001 Median :-0.253118 Median :-0.172297
## Mean : 0.003797 Mean :-0.01065 Mean : 0.006638 Mean : 0.004287
## 3rd Qu.: 0.754320 3rd Qu.: 0.82406 3rd Qu.:-0.253118 3rd Qu.:-0.172297
## Max. : 3.087256 Max. : 1.89463 Max. : 8.074589 Max. :11.368946
## alt dbarr acc
## Min. :-0.99359 Min. :-1.153377 Min. :-1.732584
## 1st Qu.:-0.61255 1st Qu.:-0.825046 1st Qu.:-0.682261
## Median :-0.33687 Median :-0.248552 Median : 0.329160
## Mean : 0.01199 Mean : 0.002765 Mean :-0.005365
## 3rd Qu.: 0.34066 3rd Qu.: 0.573661 3rd Qu.: 0.757069
## Max. : 4.41822 Max. : 3.898308 Max. : 1.846292
##
## Observation-level covariates:
## OCC
## 1 : 49951
## 2 : 49951
## 3 : 49951
## 4 : 49951
## NA's:117596
##
## Yearly-site-level covariates:
## PO year trendyear SDAC
## Min. :0.0000 01 : 3450 Min. : 1 Min. :0.000
## 1st Qu.:0.0000 02 : 3450 1st Qu.: 6 1st Qu.:0.000
## Median :0.1278 03 : 3450 Median :12 Median :0.000
## Mean :0.6243 04 : 3450 Mean :12 Mean :0.209
## 3rd Qu.:0.7983 05 : 3450 3rd Qu.:18 3rd Qu.:0.000
## Max. :8.0075 06 : 3450 Max. :23 Max. :8.000
## (Other):58650
## LDAC
## Min. :-1.0912
## 1st Qu.:-0.6831
## Median :-0.4912
## Mean : 0.0171
## 3rd Qu.: 0.5300
## Max. : 2.5113
## 3. Estimation des paramètres
Ici on ajuste le meilleur modèle d’occupancy trouvé dans le papier de Julie. Cela prend 4-5 minutes sur mon ordinateur. Un peu long, mais bcp moins que les plusieurs heures nécessaires à ajuster le modèle en bayésien.
fm <- colext(
psiformula = ~ 1, # initial occupancy
gammaformula = ~ forest + agr + halt + alt + SDAC + LDAC, # colonization
epsilonformula = ~ 1, # extinction
pformula = ~ PO + acc + OCC, # detection
umf,
control = list(trace = 1, maxit = 1e4),
se = FALSE) # pour aller plus vite, je ne calcule pas les SE et intervalles de confiance
save(fm, file = "res/dynoccloup.RData")On obtient les estimations ci-dessous.
load("res/dynoccloup.RData")
fm##
## Call:
## colext(psiformula = ~1, gammaformula = ~forest + agr + halt +
## alt + SDAC + LDAC, epsilonformula = ~1, pformula = ~PO +
## acc + OCC, data = umf, se = FALSE, control = list(trace = 1,
## maxit = 10000))
##
## Initial:
## Estimate SE z P(>|z|)
## -6.28 NA NA NA
##
## Colonization:
## Estimate SE z P(>|z|)
## (Intercept) -5.166 NA NA NA
## forest 0.929 NA NA NA
## agr 0.768 NA NA NA
## halt -0.144 NA NA NA
## alt 0.834 NA NA NA
## SDAC 0.576 NA NA NA
## LDAC 0.384 NA NA NA
##
## Extinction:
## Estimate SE z P(>|z|)
## -1.11 NA NA NA
##
## Detection:
## Estimate SE z P(>|z|)
## (Intercept) -2.497 NA NA NA
## PO 0.393 NA NA NA
## acc -0.842 NA NA NA
## OCC2 0.495 NA NA NA
## OCC3 0.436 NA NA NA
## OCC4 0.419 NA NA NA
##
## AIC: 19628.14Ces estimations sont sur l’échelle logit. Pour aller plus vite, je n’ai pas calculé les SE, d’où les NA dans les sorties du modèle. Le plus important est qu’elles sont proches des estimations obtenues en bayésien pour le papier de Julie.
4. Visualisation de l’occupancy
Pour faire une carte annuelle de la probabilité d’occupancy, soit on passe par la probabilité d’occupancy \(\psi_{i,t}\), soit on se base directement sur l’occupancy réalisée, c’est-à-dire la variable latente \(z_{i,t}\) qui nous dit si la cellule \(i\) est occupée l’année \(t\) (\(z_{i,t} = 1\)) ou pas (\(z_{i,t} = 0\)). On va se baser sur l’occupancy réalisée. Pour ce faire, il nous faut l’estimation des \(z\) qu’on obtient dans unmarked via des méthodes Bayésiennes empiriques : on estime la distribution a posteriori de la variable latente \(z\) avec les données et les estimations du max de vraisemblance obtenue au-dessus. Le mode de la distribution a posteriori est le “empirical best unbiased predictor (EBUP)”.
re <- ranef(fm)
#confint(re, level = 0.9) # 90% CI
#z_mean <- bup(re, stat = "mean") # Posterior mean
z_mode <- bup(re, stat = "mode") # Posterior mode
#head(z_mode)
#dim(z_mode)L’objet z_mode contient les \(z\) estimés avec en lignes les cellules de la grille et en colonnes les années. Mettons tout ça sur une carte. On charge une carte de la France.
pays <- st_read("shp/pays/Country.shp")## Reading layer `Country' from data source `/Users/oliviergimenez/Dropbox/Mon Mac (MacBook-Pro-de-Olivier-2.local)/Desktop/analyses/shp/pays/Country.shp' using driver `ESRI Shapefile'
## Simple feature collection with 54 features and 2 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -2177943 ymin: 5251488 xmax: 4523874 ymax: 11388590
## projected CRS: RGF93 / Lambert-93Puis notre grille complète 10km par 10km.
grid_rect <- st_read("shp/grille10par10//grille_France.shp") %>%
st_transform(crs= st_crs(pays))## Reading layer `grille_France' from data source `/Users/oliviergimenez/Dropbox/Mon Mac (MacBook-Pro-de-Olivier-2.local)/Desktop/analyses/shp/grille10par10/grille_France.shp' using driver `ESRI Shapefile'
## Simple feature collection with 5160 features and 1 field
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 510000 ymin: 6110000 xmax: 1110000 ymax: 6970000
## projected CRS: RGF93_Lambert_93On récupère les coordonnées de nos cellules échantillonnées au moins une fois au cours de l’étude.
grid <- ZST[!ind,] %>%
as_tibble() %>%
st_as_sf(coords = c('X','Y'), crs = st_crs(grid_rect))L’objet grid ainsi créé est un objet spatial sf de type POINT, il me faudrait plutôt des POLYGONS pour les représenter en couleur sur la grille.
grid_poly <- st_join(grid_rect, grid, join = st_covers)
grid_poly <- grid_poly[!is.na(grid_poly$id),]Enfin, une carte ! On prend l’année 2017 la plus récente dans le jeu de données.
grid_rect %>%
ggplot() +
geom_sf(alpha = 0, lwd = 0.01) +
geom_sf(data = grid_poly, aes(fill = as_factor(z_mode[,23])),
lwd = 0.01) +
geom_sf(data = pays %>% filter(NAME == "France"), alpha = 0) +
scale_fill_manual(name = "",
values = c("gray90",
"steelblue1"),
labels = c("cellule non-occupée",
"cellule occupée")) +
labs(title = "Carte de l'occupancy du loup en 2017",
subtitle = "Données réseau loup")
Construisons une fonction qui fait la carte de l’occupancy pour une année quelconque. Cette fonction prend comme argument year entre 1995 et 2017.
occ_plot <- function(year){
if (year < 1995) stop("Le loup venait juste d'arriver, pas la peine d'en faire une carte")
if (year > 2017) stop("Pas de données aussi récentes, va falloir attendre")
index <- year - 1994
grid_rect %>%
ggplot() +
geom_sf(alpha = 0, lwd = 0.01) +
geom_sf(data = grid_poly, aes(fill = as_factor(z_mode[,index])),
lwd = 0.01) +
geom_sf(data = pays %>% filter(NAME == "France"), alpha = 0) +
scale_fill_manual(name = "",
values = c("gray90",
"steelblue1"),
labels = c("non-occupée",
"occupée")) +
labs(subtitle = year) +
theme(legend.position = "none")
}On fait un essai avec l’année 2016.
occ_plot(2016)
On fait 3 cartes pour les années 2000, 2010 et 2017.
library(patchwork)
plot_list <- lapply(X = c(2000, 2010, 2017), FUN = occ_plot)
plot_list[[1]] | plot_list[[2]] | plot_list[[3]]
5. Comparaison ZPP et occupancy
Essayons maintenant de superposer les ZPPs sur ces cartes. On récupère les ZPPs.
zpp <- st_read("shp/ZPP2018/ZPP_HIV_2018.shp") %>%
st_transform(st_crs(pays))## Reading layer `ZPP_HIV_2018' from data source `/Users/oliviergimenez/Dropbox/Mon Mac (MacBook-Pro-de-Olivier-2.local)/Desktop/analyses/shp/ZPP2018/ZPP_HIV_2018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 112 features and 11 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 598994 ymin: 6137074 xmax: 1086276 ymax: 6868641
## projected CRS: RGF93_Lambert_93Oksana rappelle qu’il y plusieurs colonnes “X14_15” (hiver 2014), “X15_16” (hiver 2015), “ETE_2016”, “HIV_16” etc jusqu’en 2018. Pour chaque colonne il y a un code associé (chiffre entre 0 et 6) :
- Sans objet
- ZPP meute
- ZPP non meute
- Secteur à suivre (pas encore zpp)
- ZPP en sursis (1ere année sans indice si 2eme alors passe 5)
- Disparue
- ?? (idem à 5) je connais pas ce code
Donc toutes les ZPP avec un 5 dans leur colonne sont considérées comme disparue pour le suivi annuel correspondant.
Toutes les ZPPs sont-elles encore actives en 2017 ?
zpp %>%
pull(FIN_ZPP)## [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [16] NA NA NA NA NA NA NA NA NA NA NA NA 2018 NA NA
## [31] NA NA NA 2017 NA NA NA NA NA NA NA NA NA NA NA
## [46] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [61] NA NA NA NA 2018 NA NA NA 2018 2018 NA NA NA NA NA
## [76] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [91] NA NA NA NA NA NA NA NA NA 2018 NA NA NA NA NA
## [106] NA NA NA NA NA NA NAOn dirait que oui, ou bien je ne comprends pas la colonne FIN_ZPP.
Depuis quand les ZPPs sont-elles actives?
zpp %>%
select(ZPPdepuis) %>%
as_tibble() %>%
select(-geometry) %>%
count(ZPPdepuis, sort = TRUE)## # A tibble: 26 x 2
## ZPPdepuis n
## <int> <int>
## 1 NA 14
## 2 0 12
## 3 2018 11
## 4 2016 10
## 5 2017 9
## 6 2013 7
## 7 2015 7
## 8 2007 5
## 9 2005 4
## 10 2011 4
## # … with 16 more rowsJe ne comprends pas à quoi correspondents les 0 et les NA.
Bref, faisons comme si je comprenais bien ce qui est dans le jeu de données des ZPPs, on refait nos cartes avec les ZPPs.
occ_plot <- function(year){
if (year < 1995) stop("Le loup venait juste d'arriver, pas la peine d'en faire une carte")
if (year > 2017) stop("Pas de données aussi récentes, va falloir attendre")
index <- year - 1994
grid_rect %>%
ggplot() +
geom_sf(alpha = 0, lwd = 0.01) +
geom_sf(data = grid_poly,
aes(fill = as_factor(z_mode[,index])),
lwd = 0.01) +
geom_sf(data = pays %>% filter(NAME == "France"),
alpha = 0) +
scale_fill_manual(name = "",
values = c("gray95",
"steelblue1"),
labels = c("cellule non-occupée",
"cellule occupée"),
guide = guide_legend(override.aes = list(linetype = "blank",
shape = NA))) +
geom_sf(data = zpp %>% filter(ZPPdepuis <= year),
aes(color = "blue"),
show.legend = "point") +
scale_color_manual(name = "",
values = c("blue"),
labels = c("ZPP"),
guide = guide_legend(override.aes = list(linetype = "blank",
shape = 21,
size = 2))) +
labs(subtitle = year) +
theme(legend.position = "none")
}On fait un essai avec l’année 2016.
occ_plot(2016)
On fait 3 cartes pour les années 2000, 2010 et 2017.
plot_list <- lapply(X = c(2000, 2010, 2017), FUN = occ_plot)
plot_list[[1]] | plot_list[[2]] | plot_list[[3]]